function [ U, V, err, gradU, gradV ] = pmf( Y, Omega, U, V, mu, lambda)
% Omega: observed entries
% |Omega(Y - U'*V)|_F^2 + mu*|U|_F^2 + lambda*|V|_F^2

maxIter = 30000;
err = zeros(maxIter, 1);
% gradient descent
for i = 1:maxIter
    stepSize = (1e-3)/i;
    gradU = 2*V*((U'*V - Y).*Omega)' + mu*U;
    U = U - stepSize*gradU;
    
    gradV = 2*U*(Omega.*(U'*V - Y)) + lambda*V;
    V = V - stepSize*gradV;
    
    err(i) = norm((Y - U'*V).*Omega, 'fro').^2 + lambda/2*norm( U, 'fro')^2 + mu/2*norm(V, 'fro')^2; 
    
    disp(['iter' ,int2str(i),', err=',num2str(err(i))]);
    if(i > 10 && abs(err(i) - err(i-1)) < 1e-6)
        break;
    end
end
err = err(1:i);

end

